%clear all
%close all
clear all
clear global
% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the limited commitment village economy for
% both the standard specification (individual deviations) and the
% alternative specification (coalitional deviations)

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version June 2021
% =====================================================================
%====================================================
% select parameters
% ===================================================

% ======== a. individual's income process ========

village=1;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=1;                         % set to 1 for coalitional deviations
incproc=1;     Nincproc=2;                   % 1 to 4 income processes define size of measurement error in income
equalsharing=0;                    % default 0, 1 if you want to impose equal sharing in the exact solution
one_per_aut=0;                     % set this to 1 if you want deviations to imply one period of autarky
if coaldev==0
convlength=0; % the conv criterion must be fulfilled at least in convlength consecutive periods to avoid jumping 'in'
else
convlength=15;
end
filetosaveest2=['COMP_Momentsest_repl' int2str(coaldev)];
filetosave_asy=['ASYMMETRY' int2str(coaldev)];

server=0;
dataincome=1;                       % set to 1 if you want to load data moments for income, otherwise something else is done
Tessa1Tobi2=1;
buildup=0;    %%%set to 1 and change Qrho and Qbeta loop to calculate for all group sizes
if village==1
    vss_high_inddev=34-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=34-1;
elseif village==2;
    vss_high_inddev=37-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=37-1;
elseif village==3
    vss_high_inddev=30-1;                 % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=30-1;
end
if vss_high_inddev>10
    vssbuildupvec_inddev=[1:9 vss_high_inddev]; % if you choose buildup==1, rhis is the vector of village sizes for which the solution is calculated in the individual deviations case
else
    vssbuildupvec_inddev=[1:9 vss_high_inddev]
end
if server==0
    if Tessa1Tobi2==2
        outputpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/JEEA Replication/Programmes/APPROXIMATE_EXACT')
        programpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/JEEA Replication/Programmes/APPROXIMATE_EXACT')
    elseif Tessa1Tobi2==1
        
        outputpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes/APPROXIMATE_EXACT')
        programpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes')
    end
else
    outputpath=sprintf('/home/workgroups/testob')
    programpath=sprintf('/home/workgroups/testob')
end
%betavec=[0.64,0.67,0.7,0.73,0.76,0.79,0.82,0.85];
betavec=[0.86 0.88 0.9 0.92 0.94 0.96];
%betavec=[0.96]; %%just trying how high it'll go.
rhovec=1;
%main execution options
gapcritset=0.0001; %convergence criterion for change in values
gapcritsetexact=0.1;
% other execution options
Nexact=2;                         % set to 1 for approximate solution and 2 for exact
T=6200; Tinit=201;                 % number of simulation periods in total, and number of initial periods to discard
maxcerr=0.9;                          % max cons meas error in log levels
N_cerr=1;                           % number of support points for cons meas error
maxincerr=2/3;                      % max inc meas error in multiples of inc variance (in levels)
% set to one if you want to specify the rest of village in terms of average
% consumption
rov_average=1;
% set this to 1 if you want to calculate the moments for levels of
% consumption based on log consumption
momentclevinlog=1;
Rov_sb=0; % set this to 1 if you want the outside option for the Rov in the case of coalitional deviations to equal the average value from the solution to the next smaller village (second-best), rather than
% just the utility from consuming average ROV income (first-best)
% convergence criterion and maximum number of iterations
maxitgapset=400;
% this is a punishment factor: default involves a consumption penalty in the first period equal
% to caut*Pun_fac
Pun_fac=0.0;
penalty=0; % for exact solution
% grid of time-varying planner weights
%lambda=[0.1,0.125,0.15,0.175,0.2,0.25,0.3,0.35:0.05:0.65,0.7,0.75,0.8,0.825,0.85,0.875,0.9]';
lambda=[0.025:0.025:0.975]';%linspace(0.05,0.95,301)';
m=length(lambda);
if floor((m+1)/2)~=((m+1)/2)
    disp('care: lambda has to have uneven no of points in first dimension')
    stop
end
S_KEEP=nan(1000*3,2,vss_high_inddev);                  % pre-allocate the matrices for coalitional deviations, so that you don't run simulations vss times
S_rov_KEEP=nan(100,1,vss_high_inddev);              % this one is for a maximum of vss=19 for the CD
P_rov_KEEP=nan(100,1000,vss_high_inddev);
PHI_KEEP=nan(m,1000*3,2,vss_high_inddev);
IDIO_DIST_ROV_KEEP=nan(1000*3,3,vss_high_inddev);

% ====================================================
% I. Specify Income process for individual and Rest of Village (rov)
% ====================================================
% data moments for log income from ICRISAT data - all villages
% first col: Aurepalle, secdon: Kanzara, third: Shirapur

datafile='armoments_fe1';   %%%%note that after renaming, these are the unconditional income moments, not the conditional ones, but they are very similar. 

eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);

covyylag(village)=0.001;

for exact=1:Nexact
% ===================================
% iterate over income processes
% ===================================
disp('now solving for coaldev, exact - press any key to continue')
disp([coaldev exact])
%pause

Momentsest=nan(31,2,size(betavec,2),size(rhovec,2),vss_high_inddev);
indicator1=zeros(1,size(rhovec,2));
%for exact=1:3
% specify a markov process for individual income



% set minimum and maximum village size
if coaldev==1                       % for caolitional deviations,
    vss_set=[1:3]; %%%so this calculates for each group of size 2,3,4
    vss_high=vss_set(end);
else
    if buildup==0% for individual deviations, start with small village and build up, but coarser
           vss_set=[2:3]; %%%this calculates for group of size 3 and 4 (which is what we are interested in)
    vss_high=vss_set(end);
        lagind=1;
    elseif buildup==1
            vss_set=[1:3];  %%%this calculates for each group of size 2,3,4
    vss_high=vss_set(end);
        lagind=1;
        vss_high=vss_set(end);% an indicator
    end
end
Income0=nan(vss_high_inddev+1,T,size(betavec,2),size(rhovec,2),Nexact);
Csim0=nan(vss_high_inddev+1,T,size(betavec,2),size(rhovec,2),Nexact);

% ===================================
% iterate over discount factors
% ===================================
for Qbeta=1:size(betavec,2) %[1,3,5,7,9,11,13,14,2,4,6,8,10,12]%
    
    
 Qrho=1
 beta=betavec(Qbeta)
 rho=rhovec(Qrho)
        if coaldev==1
            lagind=[];
        else
            lagind=1;
        end
       
        
        P_LAG=nan(1000,1000,length(vss_set));Y_lag=nan(1000,1000,length(vss_set));
        waut_lag=[]; waut_avg_lag=[]; waut_rov_sb=[]; waut_rov=[];waut_temp=[];
         WAUT_LAG=nan(1000,1000,length(vss_set));
            STABLE=[];
        
        for vsscount=1:length(vss_set)       % iterate over village size
            %%%%assign the income values here, because you have different
            %%%%size income grid depending on whether you go up to 4 or
            %%%%just 3
            vss=vss_set(vsscount);
            vss_max=max(vss_set);
            if vsscount>1 & coaldev==1
                lagind=lagind+(vss-(vss_set(vsscount-1)+1));
            end
            
            disp('now solving for parameters')
            disp([Qbeta,Qrho])
            disp('vss')
            disp(vss)
            
            if dataincome==1
                varnu_max=maxincerr*vardy/2;                % nu is measurement error in log levels
                rho_min=covyylag(village)/vary(village);                      % AR(1) coefficient w/o meas error
                rho_max=covyylag(village)/(vary(village)-varnu_max(village));          % AR(1) coefficient w max meas error
                rho1vec=linspace(rho_min,rho_max,Nincproc);
                varnu=vary(village)-covyylag(village)./rho1vec;
                vareps=(vary(village)-varnu).*(1-rho1vec.^2);
                if  vss<=2
                [incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,3);
                elseif  vss==3
                [incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,2);
                end
                    
                
            else
                incproc=1;
                varnu=zeros(1,Nincproc);
                rho1vec=0;
                vareps=0.3;
                [incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,3);
            end




            incomevals=exp(incomevals');
            QQQ=Q1^1000;
            SS=QQQ(:,1);
            incomevals=incomevals/(SS'*incomevals);
            cumQ=cumsum(Q1,2);
           
            
            
            
            
            
            csim=[];  cc =[]; eeewaut =[]; eewaut =[]; ewaut =[]; waut =[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov =[]; P_rov =[];
                S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
                wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[]; waut_check=[]; waut_check_ind=[];
                Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
                S_rov_aut=[]; S_rov_sb=[]; S_rov_aut_all2=[]; S_rov_aut_all=[];
                w =[];  phinew =[]; Yagg =[]; wint =[]; wsb =[]; ewsb =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; ewsb =[]; bind=[];

                Nvill=ceil((vss_high_inddev+1)/(vss+1)); 
                 
            if exact==1 
                if vss<=2
                disp('now solving for approximate solution')
                IDapprox_15_Jan_20
                elseif vss==3
                P_LAG=nan(1000,1000,length(vss_set));Y_lag=nan(1000,1000,length(vss_set));
                waut_lag=[]; waut_avg_lag=[]; waut_rov_sb=[]; waut_rov=[];waut_temp=[];
                WAUT_LAG=nan(1000,1000,length(vss_set));
                STABLE=[];
                for vss=1:3
                csim=[];  cc =[]; eeewaut =[]; eewaut =[]; ewaut =[]; waut =[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov =[]; P_rov =[];
                S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
                wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[]; waut_check=[]; waut_check_ind=[];
                Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
                S_rov_aut=[]; S_rov_sb=[]; S_rov_aut_all2=[]; S_rov_aut_all=[];
                w =[];  phinew =[]; Yagg =[]; wint =[]; wsb =[]; ewsb =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; ewsb =[]; bind=[];

                Nvill=ceil((vss_high_inddev+1)/(vss+1));   
                
                 IDapprox_15_Jan_20
                end
                
                
                
                end
                if vss>1  & (STABLE(vss)==1 | coaldev==0)
                sim_approx
                end
            else
              
                if Tessa1Tobi2==2 && server==0
                cd('C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes\COALDEVEXACT')
                elseif Tessa1Tobi2==1 && server==0
                cd('C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Programmes\COALDEVEXACT')
                end
                
                disp('now solving for exact solution')
                
                vss=vss+1; 
                if vss==2
                GHMax2_15_08_2018_script
                elseif vss==3
                    if coaldev==0
                    GHMax3_Ind_dev_15_08_2018_script
                    elseif coaldev==1 
                    
                    GHMax3_Group_dev_15_08_2018_script_old
                    GHMax3_Group_dev_05_02_2020_script     
                    
                    end
                sim_exact 
                
                %%%only simulate for group size 3, don't need it for pairs
                CRITGAP(vss)=gap;
                elseif vss==4 & Qbeta>=3
                    if coaldev==0
                    GHMax4_Ind_dev_04_09_2018_script
                    elseif coaldev==1 %%%in case of 4, you need to rerun 2 and 3 again because you've changed the income grid
                    GHMax2_15_08_2018_script
                    GHMax3_Group_dev_15_08_2018_script_old
                    GHMax3_Group_dev_05_02_2020_script  
                    GHMax4_Group_dev_10_06_2021_script    
                    
                    
                    end
                                       
                    
                sim_exact
                  
                        
                CRITGAP(vss)=gap;
                end
%                
                if Tessa1Tobi2==2 && server==0
                cd('C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\Programmes')
                elseif Tessa1Tobi2==1  && server==0
                cd('C:\Users\tbold\Dropbox\Tessa and Tobi\Programmes')
                end
                vss=vss-1;
                end
            
            
            if vss==1 | (vss==3 & Qbeta<3)
            else    
                
                % =========================================================
                % calculate moments of consumption and income growth
                % =========================================================
                if coaldev==0 || exact==2 || STABLE(vss)==1
                
                csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
                rng(3); cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];logcsim=[]; yy_hc=[];ones_hc=[];yy_hc_cond=[];ones_hc_cond=[];yy_hc_groupcond=[]; ones_hc_groupcond=[];

                rng(4); incerr=randn(size(csim0));   Momentsest1=nan(27,11); ind=[];dc_log_groupcond=[];dy_log_groupcond=[];
                % loop over size of cons measurement error
               for kkk=1:size(csim0,1)
                    tempvar(kkk)=var(log(csim0(kkk,:)));
                end
                Tempvar=mean(tempvar);
                %for ic=1:N_cerr+1
                ic=1; 
                csim=exp(log(csim0)+((ic-1)/N_cerr*maxcerr/(1-(ic-1)/N_cerr*maxcerr)*Tempvar)^0.5*cerr);
                    logcsim=log(csim);
                    logcsim0=log(csim0);
                    logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
                    logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
                    % add inc meas error back in
                    income=exp(log(income0)+varnu(incproc)^0.5*incerr);
                    logincome=log(income);
                    logincome0=log(income0);
                    logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
                    income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);


                    for ivill=1:Nvill
                        aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logaggvillincome=log(aggvillincome);

                        logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);


                    end
                    aggincome=sum(income,1);



                    % === discard first Tinit-1 periods
                    % this chooses log or level consumption for
                    % moments
                    if momentclevinlog==1
                        ccc=logcsim(:,Tinit:T);
                        ccc_cond=logcsim_meandev(:,Tinit:T);
                        ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
                    else
                        ccc=(csim(:,Tinit:T));
                        ccc_cond=(csim_meandev(:,Tinit:T)); 
                    end
                    yyy=logincome(:,Tinit:T);
                    yyy_cond=logincome_meandev(:,Tinit:T);
                    yyy_groupcond=logincome_groupmeandev(:,Tinit:T);

                    cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
                    cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
                    cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);

                    yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
                    yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
                    aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));

                    
                    vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));

                    % calculate moments for each individiual
                    for ind=1:size(csim,1)

                        dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
                        dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
                        dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
                        dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));

                        % conditioning on income increases, and income
                        % decreases
                        % NB: we allocate no change to decrease!!
                        rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
                        rrhc=find(yy(ind,:)>0);
                        rrlc=find(yy(ind,:)<=0);
                        rrhc_cond=find(yy_cond(ind,:)>0);
                        rrlc_cond=find(yy_cond(ind,:)<=0);
                        rrhc_groupcond=find(yy_groupcond(ind,:)>0);
                        rrlc_groupcond=find(yy_groupcond(ind,:)<=0);


                        %%%THE VARIANCES

                        vardclog(ind,vss)=var(cc(ind,:));
                        vardylog(ind,vss)=var(yy(ind,:));
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));

                        %%%conditional on village income!
                        vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
                        vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
                        vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
                        vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));

                        %%%conditional on group income!
                        %Step 1: create your residuals
                        vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
                        vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
                        vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
                        vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));

                        %%%and the relative variances

                        reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));

                        reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
                        reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);


                        %%%THE BETA COEFFICIENTS

                        %1) Unconditional
                        %risk-sharing whole sample
                        covtemp=cov(cc(ind,:),aggyy(1,:));
                        betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
                        covtemp=cov(cc(ind,:),yy(ind,:));
                        betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY=cc(ind,:)';
                        yy_hc(ind,:)=zeros(1,size(cc,2));
                        yy_hc(ind,rrhc)=yy(ind,rrhc);
                        oneone=ones(1,size(cc,2));
                        ones_hc(ind,:)=zeros(1,size(cc,2));
                        ones_hc(ind,rrhc)=oneone(rrhc);
                        % when only two income values, income rises only
                        % take one value, so the constant multiplying
                        % ones_hc(ind,:)',yy(ind,:)' is not defined
                        if length(incomevals)>2
                        XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
                        beta_yhigh(1:4,ind,vss)=regress(YY,XX);
                        betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
                        betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
                        else
                        XX=[ones(1,size(cc,2))',yy(ind,:)',yy_hc(ind,:)'];  
                        beta_yhigh(1:3,ind,vss)=regress(YY,XX);
                        betadcdy_low(ind,vss)=beta_yhigh(2,ind,vss);
                        betadcdy_high(ind,vss)=beta_yhigh(3,ind,vss)+beta_yhigh(2,ind,vss);
                        end

                        %2) conditional on villag income
                        covtemp=[];
                        X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
                        YY=cc_cond(ind,:)';
                        beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);

                        X=[ones(size(aggyy(1,:),2),1),aggyy'];%
                        YY=cc(ind,:)';
                        beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);

                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY_cond=cc_cond(ind,:)';
                        yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
                        oneone=ones(1,size(cc_cond,2));
                        ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
                        %%%taking care of different length of incomevals
                        if length(incomevals)>2
                        XX_cond=[ones(1,size(cc,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
                        beta_yhigh_cond(1:4,ind,vss)=regress(YY_cond,XX_cond);
                        else
                        XX_cond=[ones(1,size(cc,2))',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
                        beta_yhigh_cond(1:3,ind,vss)=regress(YY_cond,XX_cond);
                        end

                        %%

                        %2) conditional on group income
                        covtemp=[];
                        X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
                        YY=cc_groupcond(ind,:)';
                        beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);

                        X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
                        YY=cc(ind,:)';
                        beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);

                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY_groupcond=cc_groupcond(ind,:)';
                        yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
                        oneone=ones(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
                        XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
                        beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);

                        %%%taking care of different length of incomevals
                        if length(incomevals)>2
                        XX_groupcond=[ones(1,size(cc,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
                        beta_yhigh_groupcond(1:4,ind,vss)=regress(YY_groupcond,XX_groupcond);
                        else
                        XX_groupcond=[ones(1,size(cc,2))',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
                        beta_yhigh_groupcond(1:3,ind,vss)=regress(YY_groupcond,XX_groupcond);
                        end







                    end
                    var_dyagglog(1,vss)=var(aggyy);
                    % vector of moments
                    % col 1-5 are parameters
                    Momentsest1(1:5,1)=[vss,beta,rho,rho1vec(incproc),varnu(incproc)]';
                    % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
                    % (8:9), conditonal on group income (10:11)
                    Momentsest1(6:11,1)=[mean(vardclog(1:size(csim,1),vss))/mean(vardylog(1:size(csim,1),vss)), mean(reldcvar(1:size(csim,1),vss))...
                        mean(vardclog_cond(1:size(csim,1),vss))/mean(vardylog_cond(1:size(csim,1),vss)), mean(reldcvar_cond(1:size(csim,1),vss))...
                        mean(vardclog_groupcond(1:size(csim,1),vss))/mean(vardylog_groupcond(1:size(csim,1),vss)), mean(reldcvar_groupcond(1:size(csim,1),vss))]';
                    % 12-14: beta, beta high, beta low unconditional
                    Momentsest1(12:14,1)=[mean(betadcdy(1:size(csim,1),vss)),mean(betadcdy_high(1:size(csim,1),vss)),mean(betadcdy_low(1:size(csim,1),vss))]';
                    % 13_16: beta, beta high, beta low conditional on
                    % village income
                    if length(incomevals)==3
                    Momentsest1(15:17,1)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(4,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)]';
                    % 18:20: beta, beta high, beta low conditional on
                    % group income
                    Momentsest1(18:20,1)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(4,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)]';
                    %21:22 are variances of dclog, dylog
                    elseif length(incomevals)==2
                    Momentsest1(15:17,1)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(2,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(3,1:size(csim,1),vss),2),mean(beta_yhigh_cond(2,1:size(csim,1),vss),2)]';
                    % 18:20: beta, beta high, beta low conditional on
                    % group income
                    Momentsest1(18:20,1)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(2,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(2,1:size(csim,1),vss),2)]';
                    %21:22 are variances of dclog, dylog
                    end
   
                        
                    Momentsest1(21:22,1)=[mean(vardclog(1:size(csim,1),vss)),mean(vardylog(1:size(csim,1),vss))];
                    %23:28, relative variance dc/dy (dy>0) and var dc/dy(dy<0);
                    %unconditinoal 23:24, conditional on village income 25:26,
                    %conditional on group income 27:28
                    Momentsest1(23:24,1)=[mean(vardc_dypos(1:size(csim,1),vss));mean( vardc_dyneg(1:size(csim,1),vss))];
                    Momentsest1(25:26,1)=[mean(vardc_dypos_cond(1:size(csim,1),vss));mean(  vardc_dyneg_cond(1:size(csim,1),vss))];
                    Momentsest1(27:28,1)=[mean(vardc_dypos_groupcond(1:size(csim,1),vss));mean( vardc_dyneg_groupcond(1:size(csim,1),vss))];

                    %%%and the critical gap
                    Momentsest1(29,1)=CRITGAP(vss);
                    %%%%and the outside option
                    Momentsest1(30,1)=mean(waut(:,1));
                    Momentsest1(31,1)=mean(waut(:,2));

                % end of ic loop - commented out end

                

                if coaldev==0
                    Momentsest(:,exact,Qbeta,Qrho,vss)=Momentsest1(:,1);
                  %  Income0(:,:,Qbeta,Qrho,vss)=income0(1:vss_high_inddev+1,:);
                  %  Csim0(:,:,Qbeta,Qrho,vss)=csim0(1:vss_high_inddev+1,:);
                elseif coaldev==1
                    Momentsest(:,exact,Qbeta,Qrho,vss)=Momentsest1(:,1);
                   % Income0(:,:,Qbeta,Qrho,vss)=income0(1:vss_high_inddev+1,:);
                   % Csim0(:,:,Qbeta,Qrho,vss)=csim0(1:vss_high_inddev+1,:);
                end
            
            
            
            
            
                end
            end
        if Qbeta==1 &  exact==2  & vsscount==2
        %%%this saves the whole workspace for beta=0.86, n=3, one version
        %%%for ID model and one version for CD model.
        cd(outputpath)
        eval(['save ', filetosave_asy, '.mat'])
        cd(programpath)   
        end    
            
        end
              


         
                toc
            
    


       Vardylog=[0.3453;0.1648;0.2786];
         Vardylog_cond=[0.336;0.162;0.268];
% 
         j=1;inddelta=6;indrho=1;vss=vss_set; 
         Momentstable(4:7,j,inddelta,indrho,vss,exact)=[Momentsest([8 15],j,inddelta,indrho,vss) ; Momentsest(9,j,inddelta,indrho,vss)/Vardylog_cond(village) ;(Momentsest([16],j,inddelta,indrho,vss)-Momentsest([17],j,inddelta,indrho,vss))];
        % pause
        filetosaveM=([filetosaveest2 int2str(exact)]);
        if buildup==0
            
           cd(outputpath)
   save(filetosaveM,'Momentsest');
cd(programpath)
        elseif buildup==1
cd(outputpath)
 save(filetosaveM,'Momentsest');
cd(programpath)
        end
    end
end


 